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Abstract 

The performance of an explicit algebraic stress model is assessed in predict- 
ing the turbulent flow and forced heat transfer in straight ducts, with square , 
rectangular , trapezoidal and triangular cross sections , under fully de veloped con- 
ditions over a range of Reynolds number's. I s€h thermal conditions are imposed 
on the duct walls , and the turbulent heat fluxes an modeled b y g radient- eliffu s i o / 1 
type models. At high Reynolds numbers (-10 5 / wall functions are used for the 
velocity and temperature 1 * * fields , while at low Reynolds numbers, damping func- 
tions are introduced into the models. Hydraulic parameters such as friction 
factor and Nusselt number are well predicted, even when damping functions are 
used , and the present formulation imposes minimal demand on the number of 
grid points without any convergence or stability problems. Comparison bet wee n 
the models is presented in terms of the hydraulic paramete rs, friction factor and 
Nusselt number , as well as in terms of the secondary flow patte rns ew cur ring 
within the ducts. 

1. Introduction 

The performance of a turbulence model in predicting t he velocity and temperature 
fields of relevant industrial problems has become increasingly important during the 

last few years. This improved predictive performance is also true for turbulent duct 
flow, which occurs frequently in many industrial applications, such as compact heat 
exchangers, gas turbine cooling systems, cooling channels in combustion chambers, 
nuclear reactors, and others. The cross section of these ducts might be both orthog- 
onal (square or rectangular) and nonorthogonal (such as trapezoidal), in which the 
generated flow is extremely complex. Sometimes, the ducts are also wavy or corru- 
gated in the streamwise direction and might be manufactured with ribs to achieve 
faster transition to turbulence. 

Several fundamental studies of turbulent flow in square and rectangular ducts exist 
in the literature. Direct numerical simulations have been carried out for a square duct 
by Gavrilakis (1992) and Huser and Biringen (1993) with Reynolds numbers of 4410 
and ~ 10 4 , respectively. Large eddy simulations for square and rectangular ducts have 
been reported by Madabhushi and Vanka (1991) at a Reynolds number of 5800, and 
by Su and Friedrich (1994) and Meyer and Rehme (1994) at Reynolds numbers up 
to 4.9 x 10 4 . Nevertheless, limitations on computational power and memory make it 
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almost impossible to directly solve for the turbulent flow field in practical engineering 
duct flows using a direct numerical simulation (DNS) approach for the foreseeable 
future. Large eddy simulations (LES) may be more tractable; although to date, their 
use has not been widespread. Thus, the prediction of the flow and heat transfer 
characteristics in engineering duct flows still requires a Reynolds-averaged approach 
using suitable turbulent closure models for both the velocity and temperature fields. 

It is known that secondary motions take place in the corner of noncircular straight 
ducts in the plane perpendicular to the streamwise flow direction. These motions are 
turbulence-induced and are commonly referred to as motions of PrandtEs second kind. 
Such motions are of importance since they redistribute the kinetic energy, influence 
the streamwise velocity, and thereby affect the wall shear stress and heat transfer. 
The effect of secondary mot ions of PrandtEs second kind on the wall shear stresses 
and heat fluxes increases considerably when the ducts are corrugated. A linear eddy 
viscosity model (LEVM) does not have the ability to predict secondary flows, but 
still it is one of the most popular models among engineers owing to its simplicity 
and overall good performance properties. Previously, Rokni and Sunden (1996, 1998) 
employed a nonlinear eddy viscosity model (NLEVM) for predicting the flow and heat 
fluxes in straight and wavy ducts with trapezoidal cross sections. This level of closure 
accounts for the Reynolds stress anisotropy arid is then able to predict the secondary 
flows within the relative cost of a two-equation formulation. 

In the study reported here, the earlier work of Rokni and Sunden (1996, 1998) 
is extended to arbitrary ducts by using an explicit algebraic stress model (EASM). 
The method is applied to square, rectangular, trapezoidal, and triangular ducts with 
iso-thermal wall conditions using gradient-diffusion type heat flux models. The heat 
flux models are a simple eddy diffusivitv model (SED), a generalized gradient dif- 
fusion hypothesis ((XIDII) model, and a model extracted from the empirical WET 
hypothesis (Launder 1988). The EASM representation is used for both low- and 
high-RevnoIds numbers without introducing any damping functions into the tensor 
representation for the Reynolds stresses; however, at low Reynolds numbers, damping 
functions are introduced into the isotropic eddy viscosity and the heat flux models. 
At high Reynolds numbers (£10 5 ), wall functions are used for both the velocity and 
temperature fields. Jayatilleke's P-function (1969) is not used because it is shown to 
be a main source of error if wall functions are used for (numerically) predicting the 
friction factor in ducts. 

One difficulty associated with turbulent convective heat transfer and fluid flow in 
ducts is obtaining satisfactory results for both friction factor and A r (/-number, if wall 
functions are used. Usually, either friction factor or A ^-number can be predicted 
satisfactorily, but not both of them. Another problem with using wall functions in 
complex geometry (duct) flows is the variability of the minimum t/ + value along the 
grid line adjacent to the boundary. (See, e.g., Rokni and Sunden, 1996.) Uniformity 
can be achieved by setting the grid points adjacent to the wall in certain positions; 
however, while this placement is easily done in orthogonal geometries, it is extremely 
difficult in nonorthogonal geometries. Alternatively, low-Reynolds number versions of 
the models usually cannot be extended to Reynolds numbers £10 4 based on hydraulic 
diameter (see e.g., Rokni and Sunden, 1999b). Therefore, it is desirable to develop a 
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model that not only predicts these hydraulic parameters satisfactorily, but that also 
can be used for a very wide range of Reynolds numbers with less demand on the 
number of grid points and without any convergence problem. 


2. Mean and Turbulent Equations 


A Reynolds-averaged Navier-Stokes (RANS) approach is used to predict the fully 
developed turbulent flow and heat transfer in the three-dimensional duct flow. The 
governing equations for the mean and turbulent quantities are 
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The turbulent stresses pr,-, (= — pv~Ti ~) and turbulent heat fluxes (pc ?) n~ 7) require 
modeling in order to close the equations. 

For the modeling of the Reynolds stresses pr,,. within the context of an algebraic 
stress formulation, transport equations for the turbulent kinetic energy and turbulent 
dissipation rate are needed: 
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where f\ = pr;jdl i/dxj is the production term. The constants (\ \ and th 2 arc* 
set to 1.44 and 1.83, respectively, and the turbulent eddy viscosity is calculated as 
Pi = pj where ~ 0.09. The functions /|, / 2 , and /,, are damping functions 
and are equal to unity in the fully turbulent region. In this study, the Abe-Kondoh- 
Nagano (1995) formulations for fi , / 2 , and /,, are used and are given by 


h = 



1 + 0.15c" 



where 



.*+ 

16 


1 + 


tie. 


o. 




\ 200 / 


( 6 ) 


Re t 



T) 


and d is the normal distance to the nearest wall. When the AI\N model is used, the 
constants <7*. and a e are both set to 1.4, and all the remaining constant coefficients are 
calibrated against the DNS channel flow data of Kim et al. (1987). At high Reynolds 
numbers, rr A . = 1 and a e is determined from j[\JC\ l {C e 2 — f 7 1 )]• 



The explicit algebraic stress model used is an extension of the Gatski and Speziale 
(1993) model and is described in Rumsey et, al. (1999). In terms of the turbulent 
stress anisotropy 6,,, it can be written as 
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where S, j and 1 1 , j are the mean strain rate and rotation rate tensors, respectively 
(S’ij + Wij = dUi/dx-j). The o,’s are scalar coefficient functions of the invariants 
if (= SijSij = {S 2 }) and ( 2 (= Bpll’,, = — {W 2 }), and are given by 
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hhe proper choice for Oj is the minimum real root of Eq. (9) (Jongen and Gatski 
1999). 

Three different models are used to provide closure for the turbulent heat flux term. 
The first is an isotropic, simple eddy diffusivitv (SED) model based on the Boussinesq 
approximation, 

//, ()T 


“jt 


o (13) 

prr T d Xj 

where the turbulent Prandtl number for temperature <tj is set to 0.89. The second is 
a model based on a generalized gradient diffusion hypothesis (GGDII), 
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and the third model is based on the WET method and is given by 
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where C, - 0.3 in both the GGDH and the WET models. At low Reynolds numbers, 
a damping function for the SED model ( J \, ) is included in the turbulent eddy viscosity 
fi t , and for the GGDH and WET models, a damping function f tl r is introduced. This 
damping function is a Lam-Bremhorst (LB) (1981) type model which is given by 

Ut = (1 - t -oo22r>Re *) 2 (1 + 2L) , (i6) 
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where Rep — py/kdjfi . If the same number of grid points was used in the cross 
section, it was found that using the LB inode] in the GGDH and WET closures gave 
better results than the AKN model for flows where Re > 10 4 . Note that the WET 
model is implicit, and the resulting system of equations for the heat fluxes are solved 
analytically in each iteration (no numerical inner iteration loop). 

Both friction factor and Nusselt number have been obtained from the computa- 
tions. The calculated friction factor is thus related to the Prandtl-law (Incropera and 
DeWitt 1996) as 



The R( number is based on the hydraulic diameter defined for two or three walls as 
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where a, 6, h, and <p are base length, upper length, height, and base angle, respectively. 
The reference to two or three walls is to the number of walls in the cross section when 
symmetry conditions are used, and .4 cross the cross-section area which can be defined 
as 0.5(<7 + h) for all cases considered here. 

The calculated A (Hiumber is related to the Ditt us- Boeder correlation (Incropera 
and DeWitt 1996) by 
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At high Reynolds numbers, the law of the wall is assumed to be valid for both 
the velocity and temperature fields in the near wall region (see Rokni and Sunden 
(1996, 1999a) for implementation details). While the log-law behavior is assumed 
for the velocity field, the temperature field can be treated by either of two methods. 
One is the usual log-law behavior, and the other is the commonly used P-function 
(Jayatilleke, 1969) method. Bv using the latter approach, the temperature field is 
given by 
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where the /’-function can be expressed as 
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and <tt = 0.89. This method is very popular, especially in commercial codes; however, 
two disadvantages of the method in duct flows are that the temperature field will he 
directly dependent on the velocity field, and that the von Karman constant k may 
play a significant roll in the determination of the N (/-number. Figure 1 clearly shows 
the effect of the k value on the calculated friction factor and Adz-number in a square 
duct (using the EASM and the SED model). A question that now arises is which 
value of the von Karman constant should be chosen. A value of k = 0.408 gives the 
best result for the friction factor; k = 0.46 gives the best result for the Ab/-number. 




Figure 1. Calculated friction factor and A r t/-number compared with experimental 
correlations in a square duct using EASM, SED, and two near wall treatments for 
t he temperature field. 

Rokni and Sunden (1996) assumed an average value of k = 0.435. In light of this 
ambiguity, the former approach of assuming a log-law behavior for the temperature 
field will be used here. 

The numerical method is based on the finite volume technique, with a nonstag- 
gered grid arrangement. The SIMPLEC algorithm is used for pressure-velocity cou- 
pling. A modified SIP method is implemented for solving the equations. The QUICK 
scheme is used for treating the convective terms in the momentum equation. How- 
ever, to achieve stability in the k and e equations, a hybrid scheme is used for the 
convective terms. A further discussion of the specification and implementation of the 
boundary conditions, as well as the numerical procedure used in the solution of the 
mean and turbulent equations, can be found in Rokni and Sunden (1996, 1999b). 

3, Results 

Straight ducts with square, rectangular, and trapezoidal cross sections, and a wavy 
duct are considered in this investigation. Only one quarter of the duct with square 
and rectangular cross sections and only half of the duct with a trapezoidal cross 
section are considered by imposing symmetry conditions. Sketches of the various 
duct configurations are shown in Fig. 2. The calculations focus on fully developed, 
three-dimensional, turbulent duct flow. Results of mean velocity, and friction factor 
and A r i/-number distributions are presented, the latter two quantities being the most 
important hydraulic parameters from an engineering standpoint. In addition, the 
secondary flow generated within the ducts is also analyzed. 


6 




NoittiWatf 


North Wall 



Figure 2. Ducts under consideration. 

3.1 Grid Sensitivity 

A nonuniform grid distribution is employed in the plane perpendicular to the main 
flow direction. Close to each wall, the number of grid points or control volumes are 
increased to enhance the resolution and accuracy. From the duct center to each wall, 
the grid distance is multiplied by a stretching factor (ST) less than unity. Thus, the 
smaller this factor, the more grid points are concentrated near the wall (i.e., finer grid 
near wall). A different number of grid points was used in the cross-sectional plane 
in order to establish the accuracy of the calculations. Table 1 shows the calculated 
A u-number and friction factor in a square duct with different stretching factors and 
number of grid points, when using the EASM. The calculated Au-number (by GGDH) 
and friction factor are compared with the correlations mentioned in previous sections. 
In the table, Nu^b stands for Au-number calculated from the Dittus-Boelter corre- 
lation, and fp r stands for friction factor calculated from the Prandtl-law correlation. 



As is evident from Table 1, decreasing the stretching factor (inserting more grid points 
in the viscous sublayers) for a specific number of grid points increases the accuracy 
of the calculations. For example, using 31 x 31 grid points with stretching factor 
0.85 yields predictions as accurate as 35 x 35 number of grid points with a stretching 
factor of 0.9. In this study, 31 x 31 grid points in the cross section with different 
stretching factors (depending on the Re number) has been used. If wall functions 


Table 1 . Calculated friction factor and N (/-number for a square duct with different 
numbers of grid points and stretching factors using low Reynolds version. 


ST 

C 

i rid 

Rt 

/ x 10 3 

J Pr X 10 3 

cliff %“ 

GGDH 

Nu db 

diff %“ 

0.9 

21 

X 

21 

4561 

10.098 

9.602 

-5.2 

16.3 

17.6 

7.4 

0.9 

:51 

X 

31 

4588 

9.981 

9.586 

-4.1 

16.2 

17.7 

8.5 

0.9 

35 

X 

35 

4603 

9.915 

9.576 

-3.5 

16.2 

17.8 

8.9 

0.9 

41 

X 

11 

4615 

9.864 

9.569 

-3.1 

16.2 

17.8 

9.0 

0.9 

51 

X 

51 

4619 

9.847 

9.567 

-2.9 

16.2 

17.8 

9.0 

0.9:} 

:$1 

X 

31 

4549 

10.151 

9.610 

-5.7 

16.6 

17.3 

4.0 

0.85 

31 

X 

31 

4604 

9.912 

9.576 

-3.5 

16.3 

1 7.8 

8.4 

0.85 

11 

X 

41 

4607 

9.900 

9.574 

-3.4 

16.3 

17.8 

8,1 


a difT ( X = 100 x (correlated - calculated )/correlated 


were used (Re > 10 1 based on hydraulic diameter), 21 x 21 to 31 x 31 (depending 
on the Rt number) grid points were sufficient to obtain reasonable accuracy for both 
the friction factor and A : (/-number in the square ducts. 

3.2 Square Duct 

The square duct is the least complicated geometry to be studied here. The flow 
and heat transfer results presented here show the wide range of Reynolds numbers 
over which the current formulation can be successfully used. 

3.2.1 Secondary Flow Patterns 

In Fig. 3, the secondary flow pattern (velocity vectors) in the fully developed re- 
gion of a square duct is shown at a Reynolds number near 4800. The results predicted 
by the EASM (with damping functions) are in excellent qualitative agreement with 
the DNS study of Monipean et al. (1996) and Gavrilakis (1992). Similar secondary 
flow patterns are predicted at both the low and high Reynolds numbers considered. 

As is well known, in laminar flow these secondary motions do not occur. In tur- 
bulent flow, the forces driving the secondary motion are concentrated in the region 
close to each corner. These motions are generated by gradients of the normal tur- 
bulent stresses. However, the linear k — e model (LEVM) (see e.g., Rokni, 1998) 
does not correctly predict these secondary motions because of its inability to accu- 
rately account for the individual normal Reynolds stresses uiul. The LEVM yields 
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Figure 3. Predicted secondary motion velocity vectors in a square duct for R( near 
4800. 

the physically incorrect expression uTi = vv = ww. It is worthwhile to point out that 
secondary motions are found with LEVM; however, they are extremely small, about 
nr 4 % to 10 3 % of the streamwise flow, and cannot normally he detected. These 
very small motions lie in the limit of numerical/computer accuracy (absolute values 
of 10“ b — 10“'). While the redistribution of the turbulent kinetic energy into the 
normal components of the Reynolds stresses is important in correctly predicting the 
secondary flow pattern, it is equally important that the Reynolds shear stress com- 
ponent be accurately predicted. The turbulent shear stress is the essential element 
in the production of the turbulent kinetic energy and as such determines the overall 
turbulent energy level of the flow. 

The predicted secondary velocity profile using the EASM, combined with both 
the wall functions (Re. — 7.1 x 10' 4 ) and damping functions (Re = 5600), is shown 
in Fig. 4. (Since the velocity vectors in the low Re case are smaller than in the 
high Re case, the results from the low Re case have been magnified (xlO) for easier 
comparison of the corresponding flow patterns.) In both Reynolds number cases, 
the secondary motions consist of two counter-rotating vortices which transport high 
momentum fluid towards the duct corner along the bisector and then outwards along 
the walls. The difference between the two Reynolds numbers lies in the spatial extent 
of the vortices within the duct. 

At low Reynolds numbers, the secondary flow close to the duct center is weak, 
and its influence on the streamwise flow is small; however, the secondary motions 
concentrated near the duct corners are strong and their effect on the streamwise flow 
is large (see Fig. 5). Nevertheless, even with the existence of the two counter-rotating 
vortices, very close to the corner, a small region of stagnant flow exists. Such a region 
is an artifact of the symmetric structure of the counter-rotating vortices in the square 
duct. As will be seen in the results from the more complicated duct geometries to be 
analyzed, such a symmetric structure no longer exists, and the corner flow pattern 
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Figure 4. Secondary motion velocity vectors predicted by the EASM with different 
near-wall treatments. 

is more complicated. Figure 5 shows the effect on the streamwise velocity contours 
[V jU bulk) predicted by the EASM using both the wall function and damping function 
approach. As can be seen for the high Reynolds number case, using wall functions 
rather than damping functions increases the predicted streamwise velocity along the 
corner bisector toward the corner. 


North Wall U-vdoaty North Wall U-veioaty 




Symmetry Wall Symmetry Wall 

Figure 5. Predicted streamwise velocity contours using EASM at two different 
Reynolds numbers and using different near-wall treatments. 


3.2.2 Hydraulic Parameters 

The accurate prediction of the friction factor and A : '(/-number is an important 
consideration in assessing turbulent model performance. In this subsection, the results 
of the computations using both wall functions and damping functions are presented. 
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The calculated friction factor and iVu-number at high Reynolds numbers using 
wall functions are shown in Fig. 6. The friction factor obtained from the EASM is 
in excellent agreement with the Prandtl-law correlation. However, the model could 
not be applied for Re numbers less than about 2.0 x 10 4 due to the large distance 
between the wall and the nearest adjacent points. 

Both the GGDH and WET methods are in excellent agreement with the Dittus- 
Boelter correlation (less than 3% deviation), while the SED method deviates some- 
what more from this correlation (sec 1 Fig. 6). The GGDH and WET models could 
not be applied for Re number less than about 2.0 x 10 4 without additional numerical 
manipulations, e.g., using the results from a higher Re number as input data for the 
lower Rt number. 




Figure 6. Calculated friction factor and A'u-number using the EASM with the wall 
functions. 


One problem associated with using the wall functions is that the grid points ad- 
jacent to a wall should be a certain distance away from the nearest wall, to get the 
average ;y + value in an acceptable range (y + > 35; see Fig. 4). The problem is more 
evident when the ducts are wavy and/or have trapezoidal cross sections. This prob- 
lem can be alleviated by using the EASM presented here. Nevertheless, the damping 
function requirement of calculating the normal distance from any point to the nearest 
wall is not an easy task in general geometries. 

Figure 7 shows that the calculated friction factor using the EASM is able to 
captures the Prandtl-law correlation (about 5% over-predicted). The figure also shows 
that the iVu-number, obtained from the GGDH and WET methods, agrees rather 
well with the Dittus-Boelter correlation, while the SED method gives less accurate 
results. It should be mentioned that the GGDH and WET models underpredict the 
Dittus-Boelter correlation for Re numbers less than about 8000. The experimental 
work of Lowdermilk et al. (1969) also shows that the Au-number in square ducts 
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Figure 7. Calculated friction factor and A r ?/-number using the EASM and damping 
functions. 

underpredicts from the Dittus-Boelter equation for Rt numbers less than about 8000. 

The presented calculation procedure is highly stable and can be extended to a 
much higher Rt number than 1CF with a minimal demand on the number of grid 
points. In Fig. 8, the calculations were performed with only 81 x 81 grid points 
for all Reynolds numbers. No convergence problems arose even at very high Rt 




Figure 8. Calculated friction factor and A r M-number using the EASM with damping 
functions at high Rt numbers. 

numbers bv using the present models. The friction factor obtained from the EASM 
is over-predicted by about 5% compared to the Prandtl-law correlation, while the 
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predicted AT/-number by the GGDH and WET closures, agrees very well with the 
Dittus-Boelter correlation. 

3.3 Rectangular Ducts 

Different rectangular ducts (side ratio 2.3,5, and 10) are considered. Even in 
these cases, if the GGDH and WET models are used, both the friction factor and 
A u-number predicted by the EASM agree very w r ell with the theoretical correlations. 
If the SED model is used, the A u-number is under-predicted by about 15%. In Fig. 
9, the secondary flow^ motion for two rectangular ducts w r ith aspect ratios of 3 and 5 
are shown. The secondary motion velocity vectors predicted by the EASM, with the 
AKN damping functions in the k — e equations, are in good agreement with what, 
has been observed in some experimental results. The existence of such secondary 
flow patterns was first, observed by Nikuradse during his experiment wit h noncircular 
ducts (see Kaka^ et ah, 1987). 


North Wall Re=11500 



0.0 0.5 1.0 1.5 2.0 2.5 3.0 

Symmetry Wall 


North Wall Re=13700 



Symmetry Wall 

Figure 9. Predicted secondary motion velocity vectors in rectangular ducts with 
aspect ratios 3 and 5 using EASM with damping functions. 


Table 2 provides calculated friction factor, AT/-number and the center-to-bulk- 
velocity ratio ( U c / 1 A ) in a rectangular duct with different aspect ratios. For a given 
cross section, the U c jUb decreases slightly with increasing Re number, which is also 
evident from this table. The experimental value of ( J c /t'b for a rectangular duct with 
aspect ratio 8 at Re ^ 5800 is 1.23 (see Rokni et ah. 1998) which can be compared 
with the calculation result (Table 2) 1.22 for aspect ratio 10 at Re « 1.572 x 10 4 . 
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Table 2. Calculated friction factor and Nu-nmnhev for rectangular ducts with differ- 
ent aspect ratios using EASM and GGDH. 


Aspect Ratio 

Re x lO- 4 

X 

O 

CJ 

Nu 

iQlh 

2 

0.9397 

8.165 

30.4 

1.28 

3 

1.1474 

7.797 

35.8 

1.28 

5 

1 .3666 

7.541 

41.2 

1.27 

10 

1.5717 

7.401 

47.3 

1.22 


3.4 Trapezoidal and Triangular Ducts 

The velocity vectors and the corresponding mean flow contours predicted by the 
EASM in a trapezoidal duct, are presented in Fig. 10. As shown in the figure, 
there exist two counter rotating vortices close to each corner that are similar to the 
results obtained by Rokni and Sunden (1996). Only 61 x 31 grid points were used 
in the cross section to perforin the calculation. Since the LB damping functions 
had convergence and stability problems regardless of grid arrangement in the cross 
section in the trapezoidal ducts, the AKN damping functions were used. In Fig. 
10. the Re number is about 1.546 x 10 4 , and the calculated friction factor and Nu- 
number (GGDH model) are 7.791 x 10~' 3 and 48. F respectively. These values can be 
compared with the Prandtl-law and Dittus-Boelter correlations, Eqs. (17) and (19), 
which yield 6.900 x 10~ 3 and 46.8. respectively. The center-to-bulk-velocity ratio 
H c/f b) is calculated as 1.29. 

Close to the upper side corner ("‘north wall" and "high wall' 1 ) there exist two 
counter-rotating vortices a small one and a much larger one. The smaller vortex 
size decreases when decreasing the upper side length ("north wall") until it vanishes 
for a triangular duct. Correspondingly, the large vortex size increases while this length 
decreases (see Fig. 11). This type of secondary flow pattern in a triangular duct was 
also observed in the experiment of Nikuradse (see e.g., Kaka^, 1987). 

The highly st able nature of the calculation procedure used here makes it possible to 
apply the present models to such triangular ducts and to predict turbulence quantities 
without any convergence problems. In Fig. 11 the upper side length is much smaller 
than the two other lengths (^ 2 x 10“ 3 of the duct height ). This length cannot be 
set to zero since using structured grids in the calculations requires that no side of 
any control volume in the domain be zero. Nevertheless, this very small upper side 
length would still yield the correct limiting behavior of a sharp corner and would be 
a case in which many turbulence models would fail. The Re number in this duct is 
1.164 x 10 4 , and the predicted friction factor and A/u-number (GGDH model) are 
8.016 x 10~ 3 and 36.3, respectively. These results can be compared with the Prandtl- 
law and Dittus-Boelter correlations, Eqs. (17) and (19), which yield 7.421 x 10 -3 and 
37.3, respectively. The center-to-bulk-velocity ratio (U c jUt>) is calculated as 1.30. 
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Re= 15500 


North Wall 




Figure 10. Velocity field in a trapezoidal duct: (a.) secondary motion velocity vectors; 
(b) mean velocity contours. 


North Wall Re=11600 North Wall Re=11600 




Figure 11. Predicted velocity and temperature fields in a triangular duct : (a) sec- 
ondary motion velocity vectors; (b)rnean temperature contours. 

3.5 Wavy Ducts 

In light of the success with the previous geometries, an initial calculation on a wavy 
duct has been done to further evaluate the performance of the model and calculation 
procedure presented in this study. The wavy duct under consideration is shown in 
Fig. 2. A symmetry plane is imposed at the cross section with an aspect ratio 4 to 3 
and sinuous variation along the //-direction. The number of grid points in the cross 
section is set to 61 x 31 for y- and c-directions, respectively. This discretization is 
similar to the number arid distribution of grid points used in the cross section for the 
trapezoidal duct. Close to each wall, the number of grid points, or control volumes, 
is increased to enhance the resolution and accuracy. Unfortunately, due to computer 
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capacity and time, only 30 grid points, uniformly spaced, are set in the streamwise, 
or .r, direction. 

For convergence, the residuals reached the value 10“ 4 for the temperature field 
and 10 -5 for the velocity field and turbulence equations. The GGDH method was 
used for the temperature equation. 

The restrictions on the streamwise resolution can adversely affect the performance 
of the solution procedure. Such inaccuracies in the computation in some regions may 
lead to large values of some key parameters which deteriorate the whole solution 
field. This situation occurs here, and to obtain a converged solution, restrictions are 
needed. The parameter 7v 2 is a useful parameter for characterizing the flow. For 
a pure shear flow Tv 2 — 1, and for a plain strain flow 1Z 2 — 0. In this study, the 
value being calculated in the cross section of the straight square duct was 0.947, for 
the trapezoidal duct, the value was 0.964, and for the triangular duct, the value was 
0.958. This range of values suggests that the model will perform well since the EASM 
was originally calibrated for homogeneous shear flows where 7Z 2 = 1. In the wavy 
wall case, values of Tv 2 exceeding 2 and correspondingly large values of 7/ greater than 
16 occur near the bend in the duct. These values yield too large values of Pk/t , which 
eventually destroy the solution. .Jongen and Gat. ski (1998) correlated the behavior of 
t hese three parameters (//.Tv 2 , Pk/t) (see their Fig. 4), arid arrived at a limiting value 
for 7v given bv 



The limiting value for Tvn, n was ±1.23. At these points Pk/t can be very large. 
Therefore, T?j- rn = 1.513 is the limiting value used in the calculations and this restricts 
the solution of the cubic equation, Eq. (9), to values of Pk/t which do not lead to a 
deteriorated solution. 

Fable 3 shows the calculated A r //-number and friction factor for the wavy duct 
in comparison with the straight trapezoidal duct. Included in the table are columns 
where the calculated friction factor has been normalized by the Prandtl-law. and the 
calculated A r //-number has been normalized by the Dittus-Boelter correlation. As can 


Table 3. Comparison between a wavy duct and straight duct with similar cross 
section. 


Type 

Rt 

/ x 10 3 

/Pr X 10 3 

./'//Fr 

GGDH 

A' 1<DB 

Nu/Nudb 

Straight. 

13699 

7.466 

7.115 

1.05 

40.4 

42.5 

0.95 

Wavy 

8944 

17.527 

7.956 

2.20 

48.2 

30.2 

1.60 


be seen from Table 3, both the friction factor and the A r //-number for the wavy duct 
are much higher than the straight duct. 

Figure 12 shows that the secondary velocity vectors at the inlet of the duct have 
changed significantly, compared to the straight duct. In addition, one can assume 
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that the secondary velocity patterns have also changed significantly in the other 
cross-section planes as well. The magnitude of these secondary flow patterns is about 
10 times larger than the secondary flow in the straight duct with similar cross section, 
or about 10% of the streamwise flow. The contours in the streamwise flow direction 


z=2 (symmetry plane) North wall 



Figure 12. (a) Secondary motion velocity vectors at the inlet of the wavy duct, 

(b) streamwise velocity contours at symmetry and middle planes, and dimensionless 
temperature contours at middle plane of wavy duct. 


are also shown in Fig. 12. The duct is moderately curved, and there is a very small 
recirculation zone in the streamwise symmetry plane of the duct near the north wall; 
however, no such recirculation exists in the middle plane. Such patterns suggest a 
complicated vortical flow field within the duct where components of vorticity in the 
cross-stream and streamwise directions may simultaneously exist. 

4. Summary 

The results from the numerical solution of fully developed, three-dimensional tur- 
bulent duct flow under isothermal conditions have been presented for square, rectan- 
gular, trapezoidal, triangular, and wavy ducts. The turbulent stresses were modeled 
using an EASM, and the turbulent heat fluxes were modeled by the SED, GGDH and 
WET methods. At high Reynolds numbers (~KF), wall functions for the velocity and 
temperature fields were used. At low Reynolds numbers, the AKN damping func- 
tions were used for the turbulent equations, and for the turbulent heat fluxes (GGDH 
and WET methods), Lam-Bremhorst type damping functions were used. Compar- 
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isons with well-established correlations, extracted from experimental studies, showed 
excellent agreement for the hydraulic parametrs (friction factor and Au-number). 
Qualitative comparisons with observed secondary flow patterns were also found to be 
in excellent agreement. 

The calculation procedure was found to be robust, with limited demand on the 
total number of grid points to achieve the desired accuracy thus minimizing the as- 
sociated computational cost. This procedure included the very challenging triangular 
duct case, where excellent results were obtained without any convergence or stability 
problems. In the wavy duct with trapezoidal cross section, streamwise resolution 
problems necessitated the imposition of a limiting value on the characteristic flow 
parameter 7 Z. Nevertheless, with this restriction, results were obtained showing the 
distinguishing features of the fully developed wavy duct flow as well as the contrasting 
behavior to the straight duct case with similar cross section. 

These results suggest that while the models for the heat fluxes can be very simple, 
this simplicity does not necessarily preclude an accurate prediction of the temperature 
field. Under isothermal conditions, simple gradient-diffusion models for the heat 
fluxes may suffice if the flow field can be well predicted. In complex geometries such 
as t hose examined here, it is necessary to use higher-order models for the Reynolds 
stresses, since anisotropies in the turbulent stress field are important, and simple 
linear eddy viscosity models will not suffice. Higher-order closures for the heat fluxes 
mav also be required in nonadiabatic cases and/or in cases where counter-gradient 
heat transfer occurs. Corresponding explicit algebraic heat flux models, coupled with 
equations for the temperature variance and variance dissipation rate, coidd be applied 
to such flows. 
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